# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# GIGER / LANZ / DEVRIES 2018

# PSRM-OA-2017-0118 -- REPLICATION MANUSCRIPT

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

rm(list=ls(all=TRUE))

library(foreign)
library(MASS)  
library(sandwich)	
library(zoo)
library(lmtest)
library(Formula)
library(betareg)
library(mfx)
library(lme4)
library(mlmRev)
library(arm)
library(stargazer)
library(RColorBrewer)
library(nloptr)
library(stringr)
library(readstata13)
library(reshape2)
library(ggplot2)
library(randomizr)

load("DATA_PSRM-OA-2017-0118.RData")

data_comp      <- data                           # complete dataset with all candidates (contacted and not contacted)
data           <- subset(data, answer!= "NA")    # data with contacted candidates (main data for analysis)
data_awn       <- subset(data, answer== "yes")   # data with contacted candidates & anwsers
data_ninc      <- subset(data, inc_r==0)         # data without incumbents
data_valais    <- subset(data, sender!=23)       # data without canton of Valais
data_staff     <- subset(data, staff==0)         # data without candidates with payed staff
data_list      <- subset(data, list_alph==0)     # data wehere listpos is non-alphabetical
data_list_alph <- subset(data, list_alph==1)     # data wehere listpos is alphabetical

options(warn=-1)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# TABLE 1: Response Rate

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# total (first variable from top)

table(data$answer)
prop.table(table(data$answer))

# incumbent (second variable from top)

table(data$answer, data$inc)
prop.table(table(data$answer, data$inc), 2)

# elected (third variable from top)

table(data$answer, data$elec)
prop.table(table(data$answer, data$elec), 2)

# gender (fourth variable from top)

table(data$answer, data$male)
prop.table(table(data$answer, data$male), 2)

#language (fifth variable from top)

table(data$answer, data$lang)
prop.table(table(data$answer, data$lang), 2)

#Political camp (sixth variable from top)

table(data$answer, data$party)
prop.table(table(data$answer, data$party), 2)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Figure 1: Kernel density distribution of answering time and length of answer

# export: 8.5in X 6in

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

par(mfrow=c(2,1))
par(mar=c(3,1.5,3,1.5), oma=c(0,0,0,0))

# Hours since contact (top panel of figure 1)

dens_h <- density(data_awn$answer_h)

plot(dens_h, ylab=" ", yaxt="n", xlim=c(0,300), ylim=c(0,0.025), 
     main="Hours since contact", type="n", 
     xaxs='i',yaxs='i', cex.main=0.85, cex.axis=0.85)

abline(v=seq(0,300,by=50), col="grey80")
abline(h=seq(0,0.03,by=0.01), col="grey80")
abline(v=seq(0,300, by=300), col="black")
abline(h=0, col="black")
points(dens_h, type="l", lwd= 2)

mean(data_awn$answer_h)
text(250, 0.022, "N=434; Mean=32.6", cex=0.85)

# Answer length (bottom chart of figure 2)

dens_l <- density(data_awn$answer_l)

plot(dens_l,  xlab="Length of answer (number of words)", 
     ylab=" ", yaxt="n",  xlim=c(0,700), ylim=c(0,0.007), 
     main="Length of answer (number of words)", type="n" , 
     xaxs='i',yaxs='i', cex.main=0.85, cex.axis=0.85)

abline(v=seq(0,700,by=100), col="grey80")
abline(h=seq(0,0.007,by=0.002), col="grey80")
abline(v=seq(0,700, by=700), col="black")
abline(h=0, col="black")
points(dens_l, type="l", lwd= 2)

mean(data_awn$answer_l)
text(600, 0.0064, "N=434; Mean=119.2", cex=0.85)

# save chart

dev.copy(pdf,"figure1.pdf", width=8.5, height=6)
dev.off()

# clean up

rm(dens_h, dens_l)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Figure 2: Treatments and response rates

# export: 8.5in X 3.75in

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

cant      <- (tapply(as.integer(data$answer), data$treatment_cant, mean)-1)*100
part      <- (tapply(as.integer(data$answer), data$treatment_party, mean)-1)*100
cant_part <- (tapply(as.integer(data$answer), data$treatment, mean)-1)*100

temp <- data.frame(
  rate  = c(cant[2], cant[1], part[2], part[1], cant_part[1], cant_part[2], cant_part[3], cant_part[4]),
  order = c(1, 2, 1, 2, 1, 1, 2, 2),
  group = c(1, 1, 2, 2, 3, 3, 3, 3),
  pos   = c(1, 1, 1, 1, 2, 3, 2, 3))

temp$group <- factor(temp$group, labels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment"))

ann_text1 <- data.frame(pos = 1, order = 1,rate = 83,lab = "In-canton",
                        group = factor("Cantonal treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))
ann_text2 <- data.frame(pos = 1, order = 2,rate = 83,lab = "Out-canton",
                        group = factor("Cantonal treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))
ann_text3 <- data.frame(pos = 1, order = 1,rate = 83,lab = "Same party",
                        group = factor("Party treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))
ann_text4 <- data.frame(pos = 1, order = 2,rate = 83,lab = "Other party",
                        group = factor("Party treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))
ann_text5 <- data.frame(pos = 2, order = 1,rate = 83,lab = "In-canton, \nsame party",
                        group = factor("Cantonal and party treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))
ann_text6 <- data.frame(pos = 3, order = 1,rate = 83,lab = "In-canton, \nother party",
                        group = factor("Cantonal and party treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))
ann_text7 <- data.frame(pos = 2, order = 2,rate = 83,lab = "Out-canton, \nsame party",
                        group = factor("Cantonal and party treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))
ann_text8 <- data.frame(pos = 3, order = 2,rate = 83,lab = "Out-canton, \nother party",
                        group = factor("Cantonal and party treatment",levels = c("Cantonal treatment", "Party treatment", "Cantonal and party treatment")))

ann_text <- rbind(ann_text1, ann_text2, ann_text3, ann_text4, ann_text5, ann_text6, ann_text7, ann_text8)

ggplot(data=temp, aes(x=order, y=rate, fill=as.factor(pos))) +
  geom_bar(stat="identity", position=position_dodge(), width=.7, colour="black") +
  facet_grid(.~ group) +
  xlab(" ") +
  ylab("Answer rate") +
  geom_text(data = ann_text, aes(label = lab), size=3, angle=90, 
            position=position_dodge(width=.9), hjust=0) +
  geom_text(aes(y=40, label=scales::percent(round(rate)/100)), 
            position=position_dodge(width=.9), 
            size=3) +
  theme(legend.position="none") +
  
  theme(panel.background = element_rect(colour = "black",fill = "white"),
        axis.line=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_text(colour="black"),
        axis.ticks=element_blank(),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.minor.y = element_blank()) +
  
  scale_y_continuous(limits=c(0,110), breaks = seq(0,80,20)) +
  scale_fill_manual(values=c("grey95", "grey75", "grey55"))

# save chart

dev.copy(pdf,"figure2.pdf", width=8.5, height=3.75)
dev.off()

# clean up

rm(temp, cant, part, cant_part,
   ann_text, ann_text1, ann_text2, ann_text3, ann_text4, ann_text5, ann_text6, ann_text7, ann_text8)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Figure 3: Cantonal treatment, electoral safety and response rates

# export: 8.5in X 3.75in

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

temp1      <- subset(data, treatment_cant== "in canton")
incanton   <- (tapply(as.integer(temp1$answer), temp1$safety_group, mean)-1)*100

temp2      <- subset(data, treatment_cant== "out canton")
outcanton   <- (tapply(as.integer(temp2$answer), temp2$safety_group, mean)-1)*100


temp3 <- data.frame(
  safety = factor(c("1","2","3","1","2","3")),  
  levels=c("1","2","3"),
  treatment = factor(c("1", "1", "1","1", "1", "1")), 
  rate = c(71.26, 72.41, 76.92, 54.48, 60.63, 63.49),
  rate = c(incanton[1], incanton[2], incanton[3], outcanton[1], outcanton[2], outcanton[3]),
  order = c(1, 1, 1, 1, 1, 1),
  group = c(1, 1, 1, 2, 2, 2))

temp3$group <- factor(temp3$group, labels = c("Treatment: in-canton", "Treatment: out-canton"))

ggplot(data=temp3, aes(x=treatment, y=rate, fill=safety)) +
  geom_bar(stat="identity", position=position_dodge(), width=.7, colour="black") +
  facet_grid(.~ group) +
  xlab(" ") +
  ylab("Answer rate") +
  geom_text(aes(y=40, label=scales::percent(round(rate)/100)), 
            position=position_dodge(width=.7), 
            size=3) +
  
  theme(panel.background = element_rect(colour = "black",fill = "white"),
        axis.line=element_blank(),
        axis.text.x=element_blank(),
        axis.text.y=element_text(colour="black"),
        axis.ticks=element_blank(),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.minor.y = element_blank()) +
  
  scale_y_continuous(limits=c(0,100), breaks = seq(0,100,20)) +
  scale_fill_manual(values=c("grey55", "grey75", "grey95"), 
                    labels=c("Low", "Medium", "High"), 
                    name="Electoral safety")

# save chart

dev.copy(pdf,"figure3.pdf", width=8.5, height=3.75)
dev.off()

# clean up

rm(temp1, temp2, temp3, incanton, outcanton)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Table 2: RI logistic regression (outcome: answer)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# m1

M1 <- glmer(answer ~ 
            +  treatment_cant_r
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M1)

# m2

M2 <- glmer(answer ~ 
            + treatment_cant_r
            + safety
            + I(treatment_cant_r*safety)
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M2)

# m3

M3 <- glmer(answer ~ 
            + treatment_party_r
            + treatment_cant_r
            + safety
            + I(treatment_cant_r*safety)
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M3)

# m4

M4 <- glmer(answer ~ 
            + treatment_party_r
            + treatment_cant_r
            + safety
            + I(treatment_cant_r*safety)
            + age
            + male_r
            + party1 + party2 + party4 + party5 
            + lang_r
            + (1 | canton), data=data, family = binomial(link = "logit"),  control=glmerControl(optimizer="nloptwrap"))
summary(M4)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# Figure 4: Responsiveness spit by cantonal treatment and across electoral safety

# export: 9in X 7in

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

par(mfrow=c(2,2))
par(mar=c(3,3,3,1.5), oma=c(1,1,0,0))

m.beta <- fixef(M4)
v.beta <- vcov(M4)

# top left panel

XX0 <- cbind(1, 
             0,
             0, 
             seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), 
             0*seq(min(data$safety[!is.na(data$safety)]), 150, by = 1),
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat0 <- XX0 %*% t(BETA)
p.hat0 <- 1/(1+exp(-y.hat0))
p.ci.XX0 <- apply(p.hat0, 1, quantile, probs=c(0.025, 0.5, 0.975))

XX1 <- cbind(1,
             0,
             1, 
             seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), 
             1*seq(min(data$safety[!is.na(data$safety)]), 150, by = 1),
             mean(data$age), 
             1,
             0,0,
             0,0,
             1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.hat1 <- XX1 %*% t(BETA)
p.hat1 <- 1/(1+exp(-y.hat1))
p.ci.XX1 <- apply(p.hat1, 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="Predicted probabilities", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="Pr(answer)", ylim=c(0.4,1), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(0.5,0.9,by=0.1), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)

text(5, p.ci.XX0[2,1]-0.02, "Out-canton", cex=0.7,  pos = 4)
text(5, p.ci.XX1[2,1]-0.02, "In-canton", cex=0.7,  pos = 4)

rug(jitter(data$safety[data$safety<=150], 0.1), ticksize = 0.02, lwd = 1)

# top right panel

fd_01 <-apply((p.hat1)-(p.hat0), 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference (in-canton vs. out-canton)", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="FD", ylim=c(-0.4,0.6), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,150, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,150,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[1,], type="l", lty="dashed", col="black")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[2,], type="l", lwd= 2, col="black")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), fd_01[3,], type="l", lty="dashed", col="black")

# bottom left panel

XX0a <- cbind(1, 
              0,
              0,
              min(data$safety[!is.na(data$safety)]), 
              0*(min(data$safety[!is.na(data$safety)])), 
              mean(data$age), 
              1,
              0,0,
              0,0,
              1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.star0a <- XX0a %*% t(BETA)
p.hat0a <- 1/(1+exp(-y.star0a))

p.hat_TEST2 <- data.frame()
for(i in 1:dim(p.hat0)[1])
{
  p.hat_TEMP <- p.hat0[i,] -  p.hat0a
  p.hat_TEST2 <- rbind(p.hat_TEST2, p.hat_TEMP)
}

p.ci.XX0 <- apply(p.hat_TEST2, 1, quantile, probs=c(0.025, 0.5, 0.975))


plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX0[2,], xlab="Electoral safety", 
      main="First difference: out-canton", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="FD",  ylim=c(-0.3,0.7), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,175, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,175,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX0[1,], type="l", lty="dashed")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX0[2,], type="l", lwd= 2)
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX0[3,], type="l", lty="dashed")

# bottom right panel

XX0b <- cbind(1, 
              0,
              1,
              min(data$safety[!is.na(data$safety)]), 
              1*(min(data$safety[!is.na(data$safety)])), 
              mean(data$age), 
              1,
              0,0,
              0,0,
              1)

BETA <- mvrnorm(1500, m.beta, v.beta)
y.star0b <- XX0b %*% t(BETA)
p.hat0b <- 1/(1+exp(-y.star0b))

p.hat_TEST2 <- data.frame()
for(i in 1:dim(p.hat1)[1])
{
  p.hat_TEMP <- p.hat1[i,] -  p.hat0b
  p.hat_TEST2 <- rbind(p.hat_TEST2, p.hat_TEMP)
}

p.ci.XX1 <- apply(p.hat_TEST2, 1, quantile, probs=c(0.025, 0.5, 0.975))

plot (seq(min(data$safety[!is.na(data$safety)]),150, by = 1), p.ci.XX1[2,], xlab="Electoral safety", 
      main="First difference: in-canton", cex.main=0.7,
      xlim=c(min(data$safety[!is.na(data$safety)]),150),
      ylab="FD",  ylim=c(-0.3,0.7), type="n", xaxt="n", yaxs="i", mgp=c(2,.5,0), cex=1, cex.axis=0.7, cex.lab=0.7)
axis(1,at=seq(0,175, by = 25), cex.axis=0.7)
abline(h=seq(-0.2,0.6,by=0.2), col="grey92")
abline(v=seq(0,175,by=25), col="grey92")
abline(h=0, col="black")

points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX1[1,], type="l", lty="dashed")
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX1[2,], type="l", lwd= 2)
points(seq(min(data$safety[!is.na(data$safety)]), 150, by = 1), p.ci.XX1[3,], type="l", lty="dashed")

# save chart

dev.copy(pdf,"figure4.pdf", width=9, height=7)
dev.off()

# clean up

rm(BETA, fd_01, p.ci.XX0, p.ci.XX1, p.hat_TEMP, p.hat_TEST2, 
   p.hat0, p.hat0a, p.hat0b, p.hat1, XX0, XX0a, XX0b, XX1, y.hat0, y.hat1, y.star0a, y.star0b,
   i, m.beta, v.beta, M1, M2, M3, M4)

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #

# END

# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #
